home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / D-G / GemsI / Original / SturmSequences.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  9.8 KB  |  580 lines  |  [TEXT/MPS ]
  1. /*
  2. Using Sturm Sequences to Bracket Real Roots of Polynomial Equations
  3. by D.G. Hook and P.R. McAree
  4. from "Graphics Gems", Academic Press, 1990
  5. */
  6.  
  7. #
  8. # Makefile
  9. #
  10. #    command file for make to compile the solver.
  11.  
  12. solve: main.o sturm.o util.o
  13.     cc -o solve main.o sturm.o util.o -lm
  14.  
  15.  
  16. /*
  17.  * solve.h
  18.  *
  19.  *    some useful constants and types.
  20.  */
  21. #define         MAX_ORDER          12    
  22. /* maximum order for a polynomial */
  23.     
  24. #define          RELERROR              1.0e-14
  25. /* smallest relative error we want */
  26.  
  27. #define          MAXPOW            32        
  28. /* max power of 10 we wish to search to */
  29.  
  30. #define          MAXIT             800        
  31. /* max number of iterations */
  32.  
  33. /* a coefficient smaller than SMALL_ENOUGH is considered to 
  34.    be zero (0.0). */
  35.  
  36. #define          SMALL_ENOUGH        1.0e-12
  37.  
  38.  
  39. /*
  40.  * structure type for representing a polynomial
  41.  */
  42. typedef      struct    p {
  43.              int    ord;
  44.              double    coef[MAX_ORDER];
  45. } poly;
  46.  
  47. extern     int        modrf();
  48. extern     int        numroots();
  49. extern     int        numchanges();
  50. extern     int        buildsturm();
  51.  
  52. extern     double    evalpoly();
  53.     
  54.  
  55.  
  56. /*
  57.  * util.c
  58.  *
  59.  *    some utlity functions for root polishing and evaluating
  60.  * polynomials.
  61.  */
  62. #include <math.h>
  63. #include <stdio.h>
  64. #include "solve.h"
  65.  
  66. /*
  67.  * evalpoly
  68.  *
  69.  *    evaluate polynomial defined in coef returning its value.
  70.  */
  71. double
  72. evalpoly (ord, coef, x)
  73.     int        ord;
  74.     double    *coef, x;
  75. {
  76.     double    *fp, f;
  77.  
  78.     fp = &coef[ord];
  79.     f = *fp;
  80.  
  81.     for (fp--; fp >= coef; fp--)
  82.     f = x * f + *fp;
  83.  
  84.     return(f);
  85. }
  86.  
  87.  
  88. /*
  89.  * modrf
  90.  *
  91.  *    uses the modified regula-falsi method to evaluate the root
  92.  * in interval [a,b] of the polynomial described in coef. The
  93.  * root is returned is returned in *val. The routine returns zero
  94.  * if it can't converge.
  95.  */
  96. int
  97. modrf(ord, coef, a, b, val)
  98.     int        ord;
  99.     double    *coef;
  100.     double    a, b, *val;
  101. {
  102.     int        its;
  103.     double    fa, fb, x, lx, fx, lfx;
  104.     double    *fp, *scoef, *ecoef;
  105.  
  106.     scoef = coef;
  107.     ecoef = &coef[ord];
  108.  
  109.     fb = fa = *ecoef;
  110.     for (fp = ecoef - 1; fp >= scoef; fp--) {
  111.         fa = a * fa + *fp;
  112.         fb = b * fb + *fp;
  113.     }
  114.  
  115.     /*
  116.      * if there is no sign difference the method won't work
  117.      */
  118.     if (fa * fb > 0.0)
  119.         return(0);
  120.  
  121.     if (fabs(fa) < RELERROR) {
  122.         *val = a;
  123.         return(1);
  124.     }
  125.  
  126.     if (fabs(fb) < RELERROR) {
  127.         *val = b;
  128.         return(1);
  129.     }
  130.  
  131.     lfx = fa;
  132.     lx = a;
  133.  
  134.  
  135.     for (its = 0; its < MAXIT; its++) {
  136.  
  137.         x = (fb * a - fa * b) / (fb - fa);
  138.  
  139.         fx = *ecoef;
  140.         for (fp = ecoef - 1; fp >= scoef; fp--)
  141.                 fx = x * fx + *fp;
  142.  
  143.         if (fabs(x) > RELERROR) {
  144.                 if (fabs(fx / x) < RELERROR) {
  145.                     *val = x;
  146.                     return(1);
  147.                 }
  148.         } else if (fabs(fx) < RELERROR) {
  149.                 *val = x;
  150.                 return(1);
  151.         }
  152.  
  153.         if ((fa * fx) < 0) {
  154.                 b = x;
  155.                 fb = fx;
  156.                 if ((lfx * fx) > 0)
  157.                     fa /= 2;
  158.         } else {
  159.                 a = x;
  160.                 fa = fx;
  161.                 if ((lfx * fx) > 0)
  162.                     fb /= 2;
  163.         }
  164.  
  165.         lx = x;
  166.         lfx = fx;
  167.     }
  168.  
  169.     fprintf(stderr, "modrf overflow %f %f %f\n", a, b, fx);
  170.  
  171.     return(0);
  172. }
  173.     
  174.  
  175.  
  176.  
  177. /*
  178.  * main.c
  179.  *
  180.  *    a sample driver program.
  181.  */
  182. #include <stdio.h>
  183. #include <math.h>
  184. #include "solve.h"
  185.  
  186. /*
  187.  * a driver program for a root solver.
  188.  */
  189. main()
  190. {
  191.     poly    sseq[MAX_ORDER];
  192.     double     min, max, roots[MAX_ORDER];
  193.     int        i, j, order, nroots, nchanges, np, atmin, atmax;
  194.  
  195.     /*
  196.      * get the details...
  197.      */
  198.  
  199.     printf("Please enter order of polynomial: ");
  200.     scanf("%d", &order);
  201.  
  202.     printf("\n");
  203.  
  204.     for (i = order; i >= 0; i--) {
  205.             printf("Please enter coefficient number %d: ", i);
  206.             scanf("%lf", &sseq[0].coef[i]);
  207.     }
  208.  
  209.     printf("\n");
  210.  
  211.     /*
  212.      * build the Sturm sequence
  213.      */
  214.     np = buildsturm(order, sseq);
  215.  
  216.     printf("Sturm sequence for:\n");
  217.  
  218.     for (i = order; i >= 0; i--)
  219.             printf("%lf ", sseq[0].coef[i]);
  220.  
  221.     printf("\n\n");
  222.  
  223.     for (i = 0; i <= np; i++) {
  224.             for (j = sseq[i].ord; j >= 0; j--)
  225.                 printf("%lf ", sseq[i].coef[j]);
  226.             printf("\n");
  227.     }
  228.  
  229.     printf("\n");
  230.  
  231.  
  232.     /* 
  233.      * get the number of real roots
  234.      */
  235.     nroots = numroots(np, sseq, &atmin, &atmax);
  236.  
  237.     if (nroots == 0) {
  238.             printf("solve: no real roots\n");
  239.             exit(0);
  240.     }
  241.  
  242.     /*
  243.      * calculate the bracket that the roots live in
  244.      */
  245.     min = -1.0;
  246.     nchanges = numchanges(np, sseq, min);
  247.     for (i = 0; nchanges != atmin && i != MAXPOW; i++) { 
  248.             min *= 10.0;
  249.             nchanges = numchanges(np, sseq, min);
  250.     }
  251.  
  252.     if (nchanges != atmin) {
  253.             printf("solve: unable to bracket all negetive roots\n");
  254.             atmin = nchanges;
  255.     }
  256.  
  257.     max = 1.0;
  258.     nchanges = numchanges(np, sseq, max);
  259.     for (i = 0; nchanges != atmax && i != MAXPOW; i++) { 
  260.             max *= 10.0;
  261.             nchanges = numchanges(np, sseq, max);
  262.     }
  263.  
  264.     if (nchanges != atmax) {
  265.             printf("solve: unable to bracket all positive roots\n");
  266.             atmax = nchanges;
  267.     }
  268.  
  269.     nroots = atmin - atmax;
  270.  
  271.     /*
  272.      * perform the bisection.
  273.      */
  274.     sbisect(np, sseq, min, max, atmin, atmax, roots);
  275.  
  276.  
  277.     /*
  278.      * write out the roots...
  279.      */
  280.     if (nroots == 1) {
  281.             printf("\n1 distinct real root at x = %f\n", roots[0]);
  282.     } else {
  283.             printf("\n%d distinct real roots for x: ", nroots);
  284.  
  285.             for (i = 0; i != nroots; i++)
  286.                 printf("%f ", roots[i]);
  287.             printf("\n");
  288.     }
  289. }
  290.     
  291.  
  292. /*
  293.  * sturm.c
  294.  *
  295.  *    the functions to build and evaluate the Sturm sequence
  296.  */
  297. #include <math.h>
  298. #include <stdio.h>
  299. #include "solve.h"
  300.  
  301. /*
  302.  * modp
  303.  *
  304.  *    calculates the modulus of u(x) / v(x) leaving it in r, it
  305.  *  returns 0 if r(x) is a constant.
  306.  *  note: this function assumes the leading coefficient of v 
  307.  *    is 1 or -1
  308.  */
  309. static int
  310. modp(u, v, r)
  311.     poly    *u, *v, *r;
  312. {
  313.     int        k, j;
  314.     double    *nr, *end, *uc;
  315.  
  316.     nr = r->coef;
  317.     end = &u->coef[u->ord];
  318.  
  319.     uc = u->coef;
  320.     while (uc <= end)
  321.             *nr++ = *uc++;
  322.  
  323.     if (v->coef[v->ord] < 0.0) {
  324.  
  325.  
  326.             for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
  327.                 r->coef[k] = -r->coef[k];
  328.  
  329.             for (k = u->ord - v->ord; k >= 0; k--)
  330.                 for (j = v->ord + k - 1; j >= k; j--)
  331.                     r->coef[j] = -r->coef[j] - r->coef[v->ord + k]
  332.                     * v->coef[j - k];
  333.     } else {
  334.             for (k = u->ord - v->ord; k >= 0; k--)
  335.                 for (j = v->ord + k - 1; j >= k; j--)
  336.                 r->coef[j] -= r->coef[v->ord + k] * v->coef[j - k];
  337.     }
  338.  
  339.     k = v->ord - 1;
  340.     while (k >= 0 && fabs(r->coef[k]) < SMALL_ENOUGH) {
  341.         r->coef[k] = 0.0;
  342.         k--;
  343.     }
  344.  
  345.     r->ord = (k < 0) ? 0 : k;
  346.  
  347.     return(r->ord);
  348. }
  349.  
  350. /*
  351.  * buildsturm
  352.  *
  353.  *    build up a sturm sequence for a polynomial in smat, returning
  354.  * the number of polynomials in the sequence
  355.  */
  356. int
  357. buildsturm(ord, sseq)
  358.     int        ord;
  359.     poly    *sseq;
  360. {
  361.     int        i;
  362.     double    f, *fp, *fc;
  363.     poly    *sp;
  364.  
  365.     sseq[0].ord = ord;
  366.     sseq[1].ord = ord - 1;
  367.  
  368.  
  369.     /*
  370.      * calculate the derivative and normalise the leading
  371.      * coefficient.
  372.      */
  373.     f = fabs(sseq[0].coef[ord] * ord);
  374.     fp = sseq[1].coef;
  375.     fc = sseq[0].coef + 1;
  376.     for (i = 1; i <= ord; i++)
  377.             *fp++ = *fc++ * i / f;
  378.  
  379.     /*
  380.      * construct the rest of the Sturm sequence
  381.      */
  382.     for (sp = sseq + 2; modp(sp - 2, sp - 1, sp); sp++) {
  383.  
  384.         /*
  385.          * reverse the sign and normalise
  386.          */
  387.         f = -fabs(sp->coef[sp->ord]);
  388.         for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
  389.                 *fp /= f;
  390.     }
  391.  
  392.     sp->coef[0] = -sp->coef[0];    /* reverse the sign */
  393.  
  394.     return(sp - sseq);
  395. }
  396.  
  397. /*
  398.  * numroots
  399.  *
  400.  *    return the number of distinct real roots of the polynomial
  401.  * described in sseq.
  402.  */
  403. int
  404. numroots(np, sseq, atneg, atpos)
  405.         int        np;
  406.         poly    *sseq;
  407.         int        *atneg, *atpos;
  408. {
  409.         int        atposinf, atneginf;
  410.         poly    *s;
  411.         double    f, lf;
  412.  
  413.         atposinf = atneginf = 0;
  414.  
  415.  
  416.     /*
  417.      * changes at positve infinity
  418.      */
  419.     lf = sseq[0].coef[sseq[0].ord];
  420.  
  421.     for (s = sseq + 1; s <= sseq + np; s++) {
  422.             f = s->coef[s->ord];
  423.             if (lf == 0.0 || lf * f < 0)
  424.                 atposinf++;
  425.         lf = f;
  426.     }
  427.  
  428.     /*
  429.      * changes at negative infinity
  430.      */
  431.     if (sseq[0].ord & 1)
  432.             lf = -sseq[0].coef[sseq[0].ord];
  433.     else
  434.             lf = sseq[0].coef[sseq[0].ord];
  435.  
  436.     for (s = sseq + 1; s <= sseq + np; s++) {
  437.             if (s->ord & 1)
  438.                 f = -s->coef[s->ord];
  439.             else
  440.                 f = s->coef[s->ord];
  441.             if (lf == 0.0 || lf * f < 0)
  442.                 atneginf++;
  443.             lf = f;
  444.     }
  445.  
  446.     *atneg = atneginf;
  447.     *atpos = atposinf;
  448.  
  449.     return(atneginf - atposinf);
  450. }
  451.  
  452. /*
  453.  * numchanges
  454.  *
  455.  *    return the number of sign changes in the Sturm sequence in
  456.  * sseq at the value a.
  457.  */
  458. int
  459. numchanges(np, sseq, a)
  460.     int        np;
  461.     poly    *sseq;
  462.     double    a;
  463.  
  464. {
  465.     int        changes;
  466.     double    f, lf;
  467.     poly    *s;
  468.  
  469.     changes = 0;
  470.  
  471.     lf = evalpoly(sseq[0].ord, sseq[0].coef, a);
  472.  
  473.     for (s = sseq + 1; s <= sseq + np; s++) {
  474.             f = evalpoly(s->ord, s->coef, a);
  475.             if (lf == 0.0 || lf * f < 0)
  476.                 changes++;
  477.             lf = f;
  478.     }
  479.  
  480.     return(changes);
  481. }
  482.  
  483. /*
  484.  * sbisect
  485.  *
  486.  *    uses a bisection based on the sturm sequence for the polynomial
  487.  * described in sseq to isolate intervals in which roots occur,
  488.  * the roots are returned in the roots array in order of magnitude.
  489.  */
  490. sbisect(np, sseq, min, max, atmin, atmax, roots)
  491.     int        np;
  492.     poly    *sseq;
  493.     double    min, max;
  494.     int        atmin, atmax;
  495.     double    *roots;
  496. {
  497.     double    mid;
  498.     int        n1, n2, its, atmid, nroot;
  499.  
  500.     if ((nroot = atmin - atmax) == 1) {
  501.  
  502.         /*
  503.          * first try a less expensive technique.
  504.          */
  505.         if (modrf(sseq->ord, sseq->coef, min, max, &roots[0]))
  506.             return;
  507.  
  508.  
  509.         /*
  510.          * if we get here we have to evaluate the root the hard
  511.          * way by using the Sturm sequence.
  512.          */
  513.         for (its = 0; its < MAXIT; its++) {
  514.                 mid = (min + max) / 2;
  515.  
  516.                 atmid = numchanges(np, sseq, mid);
  517.  
  518.                 if (fabs(mid) > RELERROR) {
  519.                     if (fabs((max - min) / mid) < RELERROR) {
  520.                         roots[0] = mid;
  521.                         return;
  522.                     }
  523.                 } else if (fabs(max - min) < RELERROR) {
  524.                     roots[0] = mid;
  525.                     return;
  526.                 }
  527.  
  528.                 if ((atmin - atmid) == 0)
  529.                     min = mid;
  530.                 else
  531.                     max = mid;
  532.             }
  533.  
  534.         if (its == MAXIT) {
  535.                 fprintf(stderr, "sbisect: overflow min %f max %f\
  536.                     diff %e nroot %d n1 %d n2 %d\n",
  537.                     min, max, max - min, nroot, n1, n2);
  538.             roots[0] = mid;
  539.         }
  540.  
  541.         return;
  542.     }
  543.  
  544.     /*
  545.      * more than one root in the interval, we have to bisect...
  546.      */
  547.     for (its = 0; its < MAXIT; its++) {
  548.  
  549.             mid = (min + max) / 2;
  550.  
  551.             atmid = numchanges(np, sseq, mid);
  552.  
  553.             n1 = atmin - atmid;
  554.             n2 = atmid - atmax;
  555.  
  556.  
  557.             if (n1 != 0 && n2 != 0) {
  558.                 sbisect(np, sseq, min, mid, atmin, atmid, roots);
  559.                 sbisect(np, sseq, mid, max, atmid, atmax, &roots[n1]);
  560.                 break;
  561.             }
  562.  
  563.             if (n1 == 0)
  564.                 min = mid;
  565.             else
  566.                 max = mid;
  567.     }
  568.  
  569.     if (its == MAXIT) {
  570.             fprintf(stderr, "sbisect: roots too close together\n");
  571.             fprintf(stderr, "sbisect: overflow min %f max %f diff %e\
  572.                 nroot %d n1 %d n2 %d\n",
  573.                 min, max, max - min, nroot, n1, n2);
  574.             for (n1 = atmax; n1 < atmin; n1++)
  575.             roots[n1 - atmax] = mid;
  576.     }
  577. }
  578.  
  579.  
  580.